Продвинутая визуализация

Task 1

Загрузите датасет very_low_birthweight.RDS
Это данные о 671 младенце с очень низкой массой тела (<1600 грамм), собранные в Duke University Medical Center доктором Майклом О’Ши c 1981 по 1987 г.
Переменными исхода являются колонки ‘dead’, а также время от рождения до смерти или выписки (выводятся из ‘birth’ и ‘exit’. 7 пациентов были выписаны до рождения).

Сделайте копию датасета, в которой удалите колонки с количеством пропусков больше 100, а затем удалите все строки с пропусками.

data <- read_rds("./very_low_birthweight.rds")
data_filtered_100 <- data[, colSums(is.na(data)) <= 100]
#data_no_na <- data_filtered_100[complete.cases(data_filtered_100),]
data_no_na <- data_filtered_100 %>% filter(rowSums(is.na(data_filtered_100)) == 0)

Task 2

Постройте графики плотности распределения для числовых переменных. Удалите выбросы, если таковые имеются. Преобразуйте категориальные переменные в факторы. Для любых двух числовых переменных раскрасьте график по переменной ‘inout’.

2.1 - Выделяем только числовые переменные и строим графики плотности распределения.

numeric_variables_only <- data_no_na[, sapply(data_no_na, is.numeric)]



for (var in names(numeric_variables_only)) {
  plot <- ggplot(numeric_variables_only, aes_string(x = var)) +
    geom_density(fill = "blue", alpha = 0.5) +
    labs(title = paste("Плотность распределения:", var), x = var) +
    theme_minimal()
  
  print(plot) 
}

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

2.2 - удаляем выбросы в количественных переменных

Используя метод IQR, я получила очень строгую фильтрацию (из 531 осталось 232 значения, включая всех dead, что важно далее. Я пробовала использовать менее строгие параметры - 2 или 3 вместо 1.5 при умножении), однако это практически не меняет ситуацию. Тут возможно не исключать строку с outlier полностью, а например заменять на NA, но я решила попробовать другой метод.

data_no_outliers <- data_no_na

for (var in names(numeric_variables_only)) {
  q1 <- quantile(data_no_outliers[[var]], 0.25, na.rm = TRUE)
  q3 <- quantile(data_no_outliers[[var]], 0.75, na.rm = TRUE)
  iqr <- q3 - q1
  upper <- q3 + 1.5 * iqr
  lower <- q1 - 1.5 * iqr
  
  data_no_outliers <- data_no_outliers %>%
    filter(data_no_outliers[[var]] >= lower & data_no_outliers[[var]] <= upper)
}

В качестве альтернативы я использовала Z-score и считала outliers только те значения, которые более чем на 3 стандартных отклонения отклоняются от среднего. Тогда из 531 значения дальше проходят 513, что кажется более разумным результатом.

data_no_outliers <- data_no_na 

for (var in names(numeric_variables_only)) {
  numeric_data <- data_no_outliers[[var]]
  z_scores <- scale(numeric_data)
  data_no_outliers <- data_no_outliers %>%
      filter(abs(z_scores) <= 3)
}

## Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
## ℹ Please use one dimensional logical vectors instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

2.3 - переводим категориальные переменные в факторы

categorical_variables <- sapply(data_no_outliers, is.character)
data_no_outliers[, categorical_variables] <- lapply(data_no_outliers[, categorical_variables], factor)

data_no_outliers %>% glimpse() %>% head()

## Rows: 513
## Columns: 19
## $ birth    <dbl> 81.514, 81.558, 81.593, 81.610, 81.624, 81.626, 81.684, 81.68…
## $ exit     <dbl> 81.539, 81.667, 81.599, 81.697, 81.700, 81.730, 81.854, 81.87…
## $ hospstay <int> 9, 40, 2, 32, 28, 38, 62, 69, 1, 93, 44, 66, 65, 44, 70, 85, …
## $ lowph    <dbl> 7.250000, 7.250000, 6.969997, 7.320000, 7.160000, 7.039997, 7…
## $ pltct    <int> 244, 182, 54, 282, 153, 229, 182, 361, 378, 255, 186, 260, 18…
## $ race     <fct> white, black, black, black, black, white, black, white, white…
## $ bwt      <int> 1370, 1480, 925, 1255, 1350, 1310, 1110, 1180, 970, 770, 1490…
## $ gest     <dbl> 32.0, 32.0, 28.0, 29.5, 34.0, 32.0, 28.0, 28.0, 28.0, 26.0, 3…
## $ inout    <fct> born at Duke, born at Duke, born at Duke, born at Duke, born …
## $ twn      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0…
## $ delivery <fct> abdominal, vaginal, abdominal, vaginal, abdominal, vaginal, v…
## $ apg1     <int> 7, 8, 5, 9, 4, 6, 6, 6, 2, 4, 8, 1, 8, 5, 9, 9, 9, 6, 1, 6, 1…
## $ vent     <int> 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1…
## $ pneumo   <int> 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0…
## $ pda      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld      <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1…
## $ year     <dbl> 81.51471, 81.55847, 81.59406, 81.61053, 81.62421, 81.62695, 8…
## $ sex      <fct> female, male, female, female, female, male, male, male, femal…
## $ dead     <int> 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0…

##    birth   exit hospstay    lowph pltct  race  bwt gest        inout twn
## 1 81.514 81.539        9 7.250000   244 white 1370 32.0 born at Duke   0
## 2 81.558 81.667       40 7.250000   182 black 1480 32.0 born at Duke   0
## 3 81.593 81.599        2 6.969997    54 black  925 28.0 born at Duke   0
## 4 81.610 81.697       32 7.320000   282 black 1255 29.5 born at Duke   0
## 5 81.624 81.700       28 7.160000   153 black 1350 34.0 born at Duke   0
## 6 81.626 81.730       38 7.039997   229 white 1310 32.0 born at Duke   0
##    delivery apg1 vent pneumo pda cld     year    sex dead
## 1 abdominal    7    0      0   0   0 81.51471 female    0
## 2   vaginal    8    0      0   0   0 81.55847   male    0
## 3 abdominal    5    1      1   0   0 81.59406 female    1
## 4   vaginal    9    0      0   0   0 81.61053 female    0
## 5 abdominal    4    0      0   0   0 81.62421 female    0
## 6   vaginal    6    1      0   0   0 81.62695   male    0

2.4 - Визуализация двух числовых переменных с раскраской по inout

ggplot(data_no_outliers, aes_string(x = "birth", y = "hospstay", color = "inout")) +
  geom_point(alpha = 0.5) +
  labs(title = "Визуализация переменных birth, hospstay с раскраской по inout") +
  theme_minimal()

В целом видим отсутствие различий в зависимости от inout.

Task 3

Проведите тест на сравнение значений колонки ‘lowph’ между группами в переменной inout. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку ‘rstatix’. Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?

Разбив значения lowph согласно переменной inout, видим, что распределения не нормальные, поэтому будем использовать непараметрический тест, в данном случае тест Манна-Уитни (Wilcoxon).

ggplot(data_no_outliers, aes(x = lowph)) +
  geom_histogram(fill = "blue", alpha = 0.5) +
  facet_wrap(~ inout) +
  theme_minimal()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

wilcox_test_out <- data_no_outliers %>%
  wilcox_test(lowph ~ inout) %>%
  add_significance()

print(wilcox_test_out)

## # A tibble: 1 × 8
##   .y.   group1       group2         n1    n2 statistic           p p.signif
##   <chr> <chr>        <chr>       <int> <int>     <dbl>       <dbl> <chr>   
## 1 lowph born at Duke transported   433    80     23662 0.000000191 ****

wilcox_test_out <- wilcox_test_out %>%
  mutate(y.position = max(data_no_outliers$lowph, na.rm = TRUE) * 1.1)

ggboxplot(data_no_outliers, x = "inout", y = "lowph", fill = "inout") +
  stat_pvalue_manual(wilcox_test_out, label = "p", tip.length = 0.01) +
  theme_minimal() +
  labs(
    title = "Сравнение групп lowph в зависимости от переменной inout"
  )

В целом, в группе transported значения ниже, чем в born at Duke. Если более низкое значение lowph ассоциировано с более низкой выживаемостью, то возможно, транспортируют более тяжелых больных, либо же сам процесс транспортировки негативно влияет на параметры пациента.

Task 4

Сделайте новый датафрейм, в котором оставьте только континуальные или ранговые данные, кроме ‘birth’, ‘year’ и ‘exit’. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.

Оставляем континуальные данные и ранговые (apg1) данные, строим корреляционную матрицу.

continuous <- data_no_outliers %>%
  select(hospstay, lowph, pltct, bwt, gest, apg1)
cor_matrix <- cor(continuous, method = "pearson", use = "complete.obs")
ggcorrplot(cor_matrix, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 3, 
           title = "Корреляции между различными переменными")

Визуализируем корреляции 2 способами - диаграмма рассеяния с линией регрессии и heatmap.

ggplot(continuous, aes(x = gest, y = bwt)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  theme_minimal() +
  labs(
    title = "Диаграмма разброса с линией регрессии для переменных с максимальной корреляцией в датасете"
    
  )

## `geom_smooth()` using formula = 'y ~ x'

heatmap(scale(continuous))

Task 5

Постройте иерархическую кластеризацию на этом датафрейме, используя euclidean distances.

continuous_scaled <- scale(continuous)
hierarchical <- hclust(dist(continuous_scaled, method = "euclidean"), method = "ward.D2")
plot(hierarchical, xlab = "", sub = "", labels = FALSE)

Task 6

Сделайте одновременный график heatmap и иерархической кластеризации. Интерпретируйте результат.

Также используем евклидово расстояние.

pheatmap(cor_matrix, clustering_distance_rows = "euclidean", 
        clustering_distance_cols = "euclidean", clustering_method = "ward.D2")

Продолжительность пребывания в госпитале обратно корррелирует с остальными показателями. В целом, это можно интерпретировать так, что чем “лучше” показатели, тем меньше требуется находиться в госпитале. Переменные bwt и gest наиболее скоррелированы. Они отражают вес при рождении и срок беременности, так что это вполне естественно.

Task 7

Проведите PCA анализ на этих данных. Проинтерпретируйте результат. Нужно ли применять шкалирование для этих данных перед проведением PCA?

pca_out <-prcomp(continuous_scaled, center = TRUE, scale. = TRUE)
summary(pca_out)

## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6
## Standard deviation     1.6049 1.0405 0.8716 0.8610 0.7590 0.51449
## Proportion of Variance 0.4293 0.1804 0.1266 0.1235 0.0960 0.04412
## Cumulative Proportion  0.4293 0.6097 0.7363 0.8599 0.9559 1.00000

summary(pca_out)$importance[2, ] 

##     PC1     PC2     PC3     PC4     PC5     PC6 
## 0.42928 0.18045 0.12661 0.12355 0.09600 0.04412

В этом случае применяем шкалирование, так как данные очень разные, имеют разную размерность. Шкалирование поможет нивелировать влияние параметров с большим абсолютным значением.

Смотрим на процент вариации, который объясняют компоненты. Первая объясняет 43%, первые три вместе - около 73%. Первые две являются наиболее важными.

Task 8

Постройте biplot график для PCA. Раскрасьте его по значению колонки ‘dead’.

biplot_graph <- cbind(continuous, pca_out$x)

biplot_graph$dead <- data_no_outliers$dead

ggplot(biplot_graph, aes(x = PC1, y = PC2, color = as.factor(dead))) +
  geom_point(alpha = 0.5) + 
  labs(title = "Biplot colored by survival",
       x = "PC1", y = "PC2") +
  scale_color_manual(values = c("blue", "red")) + 
  theme_minimal()

Не видно различий между выжившими и нет.

Task 9

Переведите последний график в ‘plotly’. При наведении на точку нужно, чтобы отображалось id пациента.

biplot_graph$id <- rownames(biplot_graph)
ggplot_biplot <- ggplot(biplot_graph, aes(x = PC1, y = PC2, color = as.factor(dead), text = id)) +
  geom_point(alpha = 0.5) + 
  labs(title = "Biplot colored by survival",
       x = "PC1", y = "PC2") +
  scale_color_manual(values = c("blue", "red")) + 
  theme_minimal()

plotly_biplot <- ggplotly(ggplot_biplot, tooltip = "text")
plotly_biplot

Task 10

Дайте содержательную интерпретацию PCA анализу. Почему использовать колонку ‘dead’ для выводов об ассоциации с выживаемостью некорректно?

Мы видим, что данные не кластеризуются, PCA не способен в данном случае различать категории выживших и нет. Некорректность выводов по колонке dead может быть связана с тем, что это бинарный показатель и он недостаточно информативен для выводов о выживаемости. Кроме того, мы только знаем, что пациент был жив в какой-то момент времени, но не знаем, что было в дальнейшем.

Task 11

Приведите ваши данные к размерности в две колонки через UMAP. Сравните результаты отображения точек между алгоритмами PCA и UMAP.

umap_result <- umap(continuous_scaled, n_neighbors = 10, min_dist = 0.1, metric = "euclidean")
biplot_graph_umap <- data.frame(UMAP1 = umap_result$layout[, 1], UMAP2 = umap_result$layout[, 2])
biplot_graph_umap$dead <- data_no_outliers$dead
biplot_graph_umap$id <- rownames(biplot_graph_umap)

ggplot_biplot_umap <- ggplot(biplot_graph_umap, aes(x = UMAP1, y = UMAP2, color = as.factor(dead), text = id)) +
  geom_point(alpha = 0.5) + 
  labs(title = "UMAP Biplot colored by survival",
       x = "UMAP1", y = "UMAP2") +
  scale_color_manual(values = c("blue", "red")) + 
  theme_minimal()

plotly_biplot_umap <- ggplotly(ggplot_biplot_umap, tooltip = "text")
plotly_biplot_umap

Как и в случае с PCA, не видим особенных различий между категориями. График выглядит слегка по-другому, имеет несколько другую форму, однако мы не видим закономерностей ни в одном, ни в другом случае.

Task 12

Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Измените основные параметры UMAP (n_neighbors и min_dist) и проанализируйте, как это влияет на результаты.

umap_result <- umap(continuous_scaled, n_neighbors = 50, min_dist = 0.1, metric = "euclidean")
biplot_graph_umap <- data.frame(UMAP1 = umap_result$layout[, 1], UMAP2 = umap_result$layout[, 2])
biplot_graph_umap$dead <- data_no_outliers$dead
biplot_graph_umap$id <- rownames(biplot_graph_umap)

ggplot_biplot_umap <- ggplot(biplot_graph_umap, aes(x = UMAP1, y = UMAP2, color = as.factor(dead), text = id)) +
  geom_point(alpha = 0.5) + 
  labs(title = "UMAP Biplot colored by survival",
       x = "UMAP1", y = "UMAP2") +
  scale_color_manual(values = c("blue", "red")) + 
  theme_minimal()

plotly_biplot_umap <- ggplotly(ggplot_biplot_umap, tooltip = "text")
plotly_biplot_umap

При изменении параметров мы видим, что график меняется. Например, теперь представители категории dead теперь скорее сгруппированы в другой его части. Также будто бы паттерн более крупный.

umap_result <- umap(continuous_scaled, n_neighbors = 5, min_dist = 0.1, metric = "euclidean")
biplot_graph_umap <- data.frame(UMAP1 = umap_result$layout[, 1], UMAP2 = umap_result$layout[, 2])
biplot_graph_umap$dead <- data_no_outliers$dead
biplot_graph_umap$id <- rownames(biplot_graph_umap)

ggplot_biplot_umap <- ggplot(biplot_graph_umap, aes(x = UMAP1, y = UMAP2, color = as.factor(dead), text = id)) +
  geom_point(alpha = 0.5) + 
  labs(title = "UMAP Biplot colored by survival",
       x = "UMAP1", y = "UMAP2") +
  scale_color_manual(values = c("blue", "red")) + 
  theme_minimal()

plotly_biplot_umap <- ggplotly(ggplot_biplot_umap, tooltip = "text")
plotly_biplot_umap

При уменьшении n_neighbors точки становятся ближе друг к другу, возможно получится различать меньшие различия.

umap_result <- umap(continuous_scaled, n_neighbors = 10, min_dist = 0.5, metric = "euclidean")
biplot_graph_umap <- data.frame(UMAP1 = umap_result$layout[, 1], UMAP2 = umap_result$layout[, 2])
biplot_graph_umap$dead <- data_no_outliers$dead
biplot_graph_umap$id <- rownames(biplot_graph_umap)

ggplot_biplot_umap <- ggplot(biplot_graph_umap, aes(x = UMAP1, y = UMAP2, color = as.factor(dead), text = id)) +
  geom_point(alpha = 0.5) + 
  labs(title = "UMAP Biplot colored by survival",
       x = "UMAP1", y = "UMAP2") +
  scale_color_manual(values = c("blue", "red")) + 
  theme_minimal()

plotly_biplot_umap <- ggplotly(ggplot_biplot_umap, tooltip = "text")
plotly_biplot_umap

При увеличении min_dist как будто бы теряется паттерн, точки более рассеянны.